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We analyze the absorption of microwaves by the Heisenberg-Ising chain combining exact calculations, based 
on the integrability of the model, with numerical calculations. Within linear response theory the absorbed 
intensity is determined by the imaginary part of the dynamical susceptibility. The moments of the normalized 
intensity can be used to define the shift of the resonance frequency induced by the interactions and the line 
width independently of the shape of the spectral line. These moments can be calculated exactly as functions of 
temperature and strength of an external magnetic field, as long as the field is directed along the symmetry axis of 
the chain. This allows us to discuss the line width and the resonance shift for a given magnetic field in the full 
range of possible anisotropy parameters. For the interpretation of these data we need a qualitative knowledge of 
the line shape which we obtain from fully numerical calculations for finite chains. Exact analytical results on the 
line shape are out of reach of current theories. From our numerical work we could extract, however, an empirical 
parameter-free model of the line shape at high temperatures which is extremely accurate over a wide range of 
anisotropy parameters and is exact at the free fermion point and at the isotropic point. Another prediction of the 
line shape can be made in the zero-temperature and zero magnetic field limit, where the sufficiently anisotropic 
model shows strong absorption. For anisotropy parameters in the massive phase we derive the exact two-spinon 
contribution to the spectral line. From the intensity sum rule it can be estimated that this contribution accounts 
for more than 80% of the spectral weight if the anisotropy parameter is moderately above its value at the isotropic 
point. 



I. INTRODUCTION 

Short-range antiferromagnetic exchange interactions are the 
predominant electronic interactions in Mott insulators. They 
are modeled by the Heisenberg-Ising Hamiltonian 
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where the sum is over nearest-neighbor sites, and J measures 
the strength of the exchange interaction. The operators s" are 
local spin- j operators, and the parameter 8 takes account of a 
possible exchange anisotropy and may include the effects of 
dipolar interactions as well. 

A sensitive experimental probe of magnetic interactions 
in solids is the absorption of microwaves, typically in ESR- 
experiments. In the simplest experimental setup a circularly 
polarized wave travels along the direction of a homogeneous 
magnetic field. The field is slowly changed and one or more 
absorption lines are observed with increasing field, whose 
precise location and width depends on the temperature. 

For experimentally accessible strengths of the incident mi- 
crowaves linear response theory 1 provides a satisfactory the- 
oretical frame for the calculation of the absorbed intensity. 
Then the key quantity to be calculated for a given Hamiltonian 
is the (imaginary part of) the dynamical susceptibility. It is 
the Fourier transform of a certain dynamical spin-spin corre- 
lation function (see below). Since such a quantity cannot be 
calculated for an interacting many-body system as the anti- 
ferromagnetic Heisenberg-Ising model <TTT>, various kinds of 
approximations have been tried in the past. Most of these ap- 
proximations break down, when many-body correlation effects 
are strong, in particular in one- and two-dimensional systems 
with strong exchange interactions. 



In this work we shall concentrate on the one-dimensional 
spin- j case which is not covered by the more traditional 
approaches. 1 2 It is relevant for the description of quasi one- 
dimensional compounds^ with strong exchange interactions. 
This case has been successfully studied by field theoretical 
methods 5 which are, however, restricted to small temperatures 
and to a limited range of magnetic fields. They also seem to 
have built in certain a priori assumptions about the line shapes. 
In one dimension purely numerical approaches 6 7 are efficient 
as well. They are unbiased, yet the extrapolation of the data to 
the thermodynamic limit of large chains may require additional 
justification and support from analytical calculations. 

The aim of the present work is to establish a number of exact 
results for the microwave absorption of the one-dimensional 
spin- j Heisenberg-Ising magnet, with the homogeneous mag- 
netic field along the magnetic symmetry axis of the chain, and 
to interpret these results in the light of numerical calculations. 
In turns the quality and validity of the numerical calculations 
can be estimated from the analytical results. 

Our work is motivated by the recent progress in calculating 
static short-range correlation functions of the integrable spin- j 
Heisenberg-Ising chain at finite temperatures and magnetic 
fields. This makes it possible to extend a remarkable result 
for the resonance shift in one-dimensional antiferromagnetic 
chains, which was obtained by Maeda et al. 8 and which utilizes 
the exact nearest-neighbor correlation functions of the isotropic 
spin- j Heisenberg chain. We present an alternative framework 
for the derivation of the resonance shift which, in the limit 
of small anisotropy, reproduces the result of Ref. 8 In our 
approach the anisotropy is treated non-perturbatively. It allows 
us, moreover, to derive an exact formula for the line width at 
fixed magnetic field. 

We utilize the fact that the absorbed intensity I(co,h) is a 
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positive function, whose integral over CO exists. The field- 
dependent moments of the corresponding normalized intensity 
turn out to be static short-range correlation functions which can 
be calculated directly for the infinite chain. The first moment 
is the average absorption frequency. In case that there is a 
single pronounced absorption line it gives a measure for the 
shift of the resonance compared to the paramagnetic absorp- 
tion frequency co = h (in the units used in this work). This 
measure is completely independent of the actual shape of the 
spectral line. Similarly, the second moment provides a shape- 
independent measure of the line width. The idea of using 
moments was introduced by van VlecP even before the linear 
response theory was created. Here we combine it with the 
finite-temperature linear response theory. The results of van 
Vleck are then recovered in the infinite-temperature limit. 

Alternatively we may normalize the intensity by its inte- 
gral over h. In this case the moments cannot be expressed 
by finite-range correlation functions. Still, these frequency- 
dependent moments can be expanded into an infinite series of 
field-dependent moments®! which is a useful starting point for 
approximations such as the high-temperature expansion or an 
expansion for small anisotropy. Interestingly, the line width 
defined in terms of the frequency-dependent moments shows a 
temperature behavior rather different from that determined by 
the field-dependent moments. 

When there is more than a single resonance, the interpreta- 
tion of the moments is less intuitive. In case of two peaks, for 
instance, the first moment would be something like the average 
location of the two peaks. For this reason it is desirable to 
have some knowledge about the full absorption spectrum ('the 
line shape'). Hence, we complemented our exact calculation 
of the moments with numerical calculations of the dynamical 
susceptibility on finite lattices up to 32 sites. The combina- 
tion of both, the exact calculation of the moments and the 
numerics, allows us to propose a model for the line shape in 
the high-temperature limit which has no free parameters. The 
actual parameters of the corresponding distribution function 
(a normal-inverse Gaussian) are determined from the first four 
exactly calculated moments. 

An unbiased but approximate calculation of the line shape 
is possible in the massive ground state phase of the model 
at vanishing magnetic field. For an isotropic system there is 
no absorption without an external field. For sufficiently high 
anisotropy, however, the absorption becomes large. In the 
massive phase the matrix elements of the local spin operators 
between the ground states and excited states ('form factors') 
are exactly known 1 1 and generally non-vanishing in the thermo- 
dynamic limit. They are classified as 2n-spinon contributions 
according to the (even) number of elementary excitations in- 
volved. Here we calculate the two-spinon contribution to the 
absorbed intensity exactly. From the intensity sum rule we 
infer that for anisotropics moderately above the isotropic point 
the two-spin contribution is dominant and amounts to more 
than 80% of the absorbed intensity. For growing anisotropy it 
rapidly approaches 100%. But as opposed to the situation with 
the dynamic structure factor for which the relative contribution 
of the two-spinon excitation is still dominant in the isotropic 
limitj 12 " 14 4t drops off rapidly for the dynamical susceptibility. 



II. THE METHOD OF MOMENTS 

For any spin system with Hamiltonian H linear response 
theory relates the intensity absorbed from a circularly polarized 
electro-magnetic wave, whose wave length is large compared 
to the distance between the spins, to the (imaginary part of the) 
dynamical susceptibility^ 

X'i-{co,h) = ^j°Jh#»{[S±(t),!T]) T . (2) 

Here S ± = S x ± iS\ and the S a = £^ =1 sf, a = x,y,z, are 
the components of the total spin. L is the number of lat- 
tice sites, (-)t stands for the canonical average at tempera- 
ture T calculated by means of the statistical operator p = 
^-(H-hs^lT i Xx ^-(H-hs z )/T^_ Through this average the dy- 
namical susceptibility depends on T and on an external ho- 
mogeneous magnetic field h which is usually applied in ESR 
experiments. The time evolution of S + in |2]) must be calcu- 
lated with H — hS z . The absorbed intensity per spin, normal- 
ized by the intensity of the incident wave and averaged over a 
half -period njco of the microwave field, is 

I(co,h) = ® X 'l-(co,h). (3) 

In order to keep this paper self-contained we included a deriva- 
tion of Q and of Q in App. [A] 

In this work we shall exclusively concentrate on the one- 
dimensional version of the Heisenberg-Ising (orXXZ) Hamilto- 
nian ([TJ. This Hamiltonian is in the class of integrable Hamilto- 
nians, but so far this does not mean that dynamical correlation 
function such as the expectation value ([S + (f),S~]) r in |2]) 
could be calculated analytically. We are only aware of three 
very special cases where this is possible. These are the 'free 
Fermion case' 8 = — 1 at T —> °°, the isotropic limit 5 — s- 
and the Ising limit 8 — > °°. We shall discuss these cases be- 
low. In the general case so far only the short time behavior 
of ([^(f),^ - ])^^^ can be accessed by directly calculating 
the commutators involved in the time evolution up to a certain 
power. We generated this series up to the order f 38 (cf. App.[c|. 
Still, the results for the short-time behavior alone are not help- 
ful for calculating the right hand side of 

Interestingly enough, as we have shown, 10 some more el- 
ementary spectral characteristics, such as the position of the 
resonance or the line width, are easier to calculate. They may 
be expressed in terms of certain static correlation functions 
that determine the moments of a normalized intensity function 
in one variable. Since I(co,h) is a function of the frequency co 
and of the magnetic field h we may normalize it by dividing 
either by the integral over CO or by the integral over h. 

In the first case we interpret the resulting normalized func- 
tion as a distribution function of frequencies which depends 
parametrically on the magnetic field. Then its moments m n are 
field dependent. This corresponds to an experimental situation 
where the field is kept fixed and the frequency is varied. We 
shall see that, from a theoretical perspective, this case is com- 
paratively simple, since the moments depend only on static 
correlation functions of finite range. The lowest moments 
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which determine the position of the resonance and its width 
can be expressed by correlation functions ranging over up to 
four lattice sites, which can be calculated exactlypSEl 

In the second case, when the intensity is normalized as a 
function of the magnetic field, the corresponding moments M„ 
depend on the frequency. This corresponds to the standard ESR 
setup in which the magnetic field is slowly swept at fixed fre- 
quency. As we shall see below this case is more sophisticated 
for a theoretical analysis, since static correlation functions for 
arbitrary distances are already involved in the calculation of the 
lowest moments. Still, the M„ can be expanded into an infinite 
series in terms of the moments m„ and their derivatives, which 
may serve as a starting point for systematic approximations. 



A. Field-dependent moments 

We temporarily assume that our chain is large but finite. 
Then the spectrum is bounded and the integrals 



The first few of them can be easily calculated. In particular, 
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exist for all non-negative integers n. Since I(co,h) is non- 
negative everywhere and since Iq > 0, we may interpret 
I((0,h)/Io as a probability distribution with moments /„. As 
we shall see, in our case it is convenient to express the /„ in 
terms of another closely related sequence of integrals 



m n (T,h) =J 



da> 

,2k 



(co-h) n X 'l-{(o,h) 
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which, by slight abuse of language, we shall call (shifted) 
moments as well. The existence of the integrals is obvious for 
every finite chain. 

Now by definition the shift of the resonance for fixed h is 
the deviation of the average frequency from the paramagnetic 
resonance at co = h, 



I\ Jni2 + hmi 

8co = h =J — 

Iq J mi +/imo 



(6) 



A measure for the line width is the mean square deviation from 
the average frequency 



h 
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A ® = T- - rpi = J' T 

Iq /q Jmi+hmo 



-8co 2 
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We see that, in order to calculate the resonance shift and the 
line width, we need to know the first four shifted moments mo, 
mi, rri2, of the dynamic susceptibility %+-■ 

In the following we shall employ the notation adx • = [X,-] 
for the adjoint action of an operator X. Then S + (t) = 
e- iht e ,adH S + , since [H,S Z ] = and [S Z ,S+] = S + , and it fol- 
lows with Q and Q that 
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(8) 



The latter formula shows that the moments m n are static 
correlation functions whose complexity grows with growing n. 



(9) 



which is the magnetization per lattice site. For the subsequent 
moments, which do not have such an immediate and simple 
interpretation, we obtain 



mi = 8(s^s 2 — 2s\s2)t 
2 

mj, = -8 2 (2s^S2S^S^ +45^2 S 3 S 4 ~^ s i s z s 3 s 4 



ni2 = -S (,Sj +4444 — 4*1*2 i 3 )t ■ 



(10a) 
(10b) 



+ 5(8^4^4 + 25+^ - 844)> r . (10c) 

These moments are certain combinations of static short-range 
correlation functions which implies, in particular, that they all 
exist in the thermodynamic limit L — >• °°. Hence, we may relax 
our restriction that we are dealing with a finite chain at this 
point. An interesting conclusion which can be drawn from 
the existence of the moments is that the field-dependent line 
shape cannot be Lorentzian, as is sometimes assumed in the 
literature, since the second moment of a Lorentzian does not 
exist. In fact, in our numerical data for finite chains we see an 
exponential decay away from the resonance (see below). Note 
that for finite magnetic field mo is of order 1, mi is of order 8, 
n%2 is of order S 2 , but all higher moments are of order S 2 as 
well. This is clear from ([8]) and will be relevant below. 

The formulae ( fT0| ) are appealing from a theoretical per- 
spective, since, due to recent progress in the theory of in- 
tegrable systems, static short-range correlation functions of 
the Heisenberg-Ising chain can be calculated exactly. It has 
been shown that all static correlation functions of the one- 
dimensional Heisenberg-Ising model are polynomials in the 
derivatives of three functions 17 <p, co, and co' which, as is com- 
mon in integrable models, can be expressed in terms of the 
solutions of certain numerically well behaved linear and non- 
linear integral equations P^l We provide the definition of these 
functions in the critical case (-1 < 8 < 0) in App.|D] For the 
massive case (8 > 0) the definitions are similar and can be 
found in Ref.[T6l 

Although, in principle, all static correlation functions of 
the Heisenberg-Ising chain in the thermodynamic limit are 
known exactly, their explicit calculation works out only at short 
distances. At larger distances the amount of computer algebra 
involved in the calculations grows excessively. In Refs. [l9l[T5l 
and [16] we obtained all correlation functions ranging over at 
most four lattice sites. This is just enough to calculate the 
moments mo,mi,m2,m3. For the simpler case of the isotropic 
model in vanishing magnetic field we obtained the correlation 
functions ranging over up to seven lattice sitesp^l 

When considering the Heisenberg-Ising Hamiltonian as an 
integrable model it is customary to parameterize all func- 
tions by a deformation parameter q in terms of which the 
anisotropy is 8 — (q — l) 2 /2q. Employing the shorthand 
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notations <p (n) = d£<p{x) \ x =o, /(„,,„) = d™d?f(x,y) \ x=y=0 , for f = CO, CO 1 , we obtairP 
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+ 1 92q 4 {q 2 -l) 2 (q 4 -l)(q 6 + 18q 4 + 8q 3 + 18q 2 + l) CO {0j2) 

+ 64q 2 (q 2 - 1 ) 2 (q 4 - 1 ) (q + 1 ) 2 (2 9 8 - 5? 7 + 26? 6 - 49q 5 + 28q 4 - 49q 3 + 26q 2 - 5q + 2)fi> (1)1) 

- 16<? 2 (<? 2 - 1 ) 2 (<? 4 - 1 ) (q + 1 ) 2 (<? 8 - q 1 + q 6 + q 5 + 2q 4 + q 3 + q 2 - q + 1 ) (2(» ( , . 3) - 3© (2 , 2) ) 
+ 64^ 2 (q 4 - 1 ) ( ? 6 - 1 ) (3? 8 + 2q 6 + 24q 5 - 1 30q 4 + 24q 3 +2q 2 +3) co' {01) 

+ (q 4 -l){q 6 -l){q+l) 2 {q 1Q -2q 9 + 25q* + I6q 7 + 1 1 8q 6 + I64q 5 

+ 1 18q 4 + I6q 3 + 25q 2 - 2q + 1 ) (<o[ 3) 0)' {l 2) + Co[ Q 1} £0 ( ' 2j3) ) 

- 1 536<? 5 (<? 4 - 1 ) (4<? 8 - 9<? 7 - 2q 6 -6q 5 +8q 4 - 6q 3 -2q 2 -9q + 4) £B (0)0) 

+ 4 ? 2 ( ? 6 - 1 ) [q + 1 ) 2 {q 2 + 1 ) (5? 8 - 2q 1 + 32q 6 + 50q 5 + 70q 4 + 50q 3 + 32q 2 -2q + 5) 

( 2 »(1,3) W (0,l) ~ 3«(2,2)®(0,1) + m {0,2)<o[ Qt3 ) - 2a) (ljl )£»[ 0|3) ~ 3(8(0,2) t»(i )2 ) " ®(0,0) ®(2,3) ) 

- 16g 2 (tf + 1 ) 2 (? 16 - ? 15 + 8^ 14 + 9^ 13 + 47 ? 12 + 45q U +96q l0 +9lq 9 + I28q & +9lq 7 + 96q 6 
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These are the moments in the thermodynamic limit. We can 
calculate them with high numerical accuracy over the whole 
range of the phase diagram, for all temperatures and magnetic 
fields as well as for arbitrary anisotropy 5. In Figs. TJ3 we 
show two examples not too far away from the isotropic point, 
namely 8 = —0.1 in the critical phase and 5 = 0.25 in the 
massive phase. Values not too far away from the isotropic point 
are most relevant for real materials. In both cases we observe 
an increase of the resonance shift 8co and a broadening of the 
spectral lines, measured as an increase of A©, for decreasing 
temperatures. 

This is interesting as it seems to contradict experimental 
results^ which claim a narrowing. This discrepancy is due to 
the different measures for the line width here and in the experi- 
mental literature. The mean square deviation used in Figs.[T]|3] 
is a customary measure for the width of wave functions in quan- 
tum mechanics. For Gaussians it is of the order of magnitude 
of an intuitive line width drawn by eye, but for distributions 
which have long and shallow tails this is no longer the case. 
Hence, Aco may strongly deviate from a typical measure of the 
line width used in the interpretation of experimental data as 



e.g. the distance between the inflection points right and left 
to the maximum of the intensity ('peak-to-peak width'). This 
discrepancy was already noted by van VleckPOne advantage 
of the mean square deviation from the resonance frequency as 
a measure of the line width is that it is defined independently of 
the line shape. In principle it should be no problem to extract it 
from experimental data. Yet, we expect that in cases, where the 
contributions from the tails of the spectral line are important, a 
problem might be to resolve these tails from the 'background' . 

The difference between different measures of the line width 
becomes rather clear from our numerical analysis below (com- 
pare also Ref. [7). For high temperatures, where we could 
extract a model for the line width from our numerical data, 
the 'peak-to-peak width' is much smaller than Aco. This can 
be attributed to the shallow tails of the absorbed intensity. In 
experiments these tails may be misinterpreted as stemming 
from couplings of the spin chain to other degrees of freedom 
and may lead to an overestimation of the background. On the 
other hand, tails are expected to have less influence on the 
resonance shift. As long as they are not too asymmetric the 
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FIG. 1. Resonance shift 8co/J and line width Aco/J in the critical 
regime at 5 = —0. 1 as function of the magnetic field. Crosses from 
fully numerical calculation for finite chain Hamiltonians of 16 and 24 
sites. 
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FIG. 2. Resonance shift 8(0 /J and line width Aco/J in the critical 
regime at 8 = —0. 1 as function of the temperature. Crosses from fully 
numerical calculation for finite chain Hamiltonians of 16 and 24 sites. 



shift 8 co of the average of the absorbed intensity should agree 
with the shift of its maximum, which is the common measure 



FIG. 3. Resonance shift 8co/J and line width Aco/J in the massive 
regime at 8 = 0.25 as function of the temperature. 



in experiments. 

For other values of the anisotropy parameters the resonance 
shift and the line width in the critical regime for, — 1 < 8 < 0, 
show a qualitatively similar behavior as in Figs. [T] and [2] Both, 
8 CO and A©, increase with decreasing temperature at fixed 
magnetic field. Note that in the massive regime, exemplified 
with Fig. [3] the resonance shift may behave non-monotonically 
as a function of temperature. 

B. Frequency-dependent moments 

In order to obtain the resonance shift and the line width 
as defined in the previous section experimentally one would 
have to measure the microwave absorption at fixed Zeeman 
field h for various values of the frequency CO and then calculate 
the required averages as integrals over co. In current ESR 
experiments different data sets are recorded. The microwave 
frequency co is kept fixed and the absorbed intensity I(co,h) = 
COx1-((0,h)/2 is determined as a function of h. This intensity 
function can be normalized by dividing by its /i-integral, and 
the corresponding frequency-depending moments define the 
resonance shift and line width in '/i-direction'. Away from 
the isotropic point (8 = 0), where /i) is symmetric and 

the absorption line is extremely narrow, these may be rather 
unrelated to their field-dependent counterparts of the previous 
section. 

In analogy with <|5j we define the frequency-dependent mo- 
ments 

M„ (T, CO) = J'" f^(h- (Ofx'i- • d2) 
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These can be expressed in terms of the m„ and their derivatives. 
Denoting the A:th derivative with respect to the second argument 
by a superscript (k) we obtain the representation 



m„(t» = (-i; 



k=0 



{-J)\ 
k\ 



„(T,(0), 



(13) 



for the frequency-dependent moments. This representation 
involves static correlation functions for arbitrarily large dis- 
tances. For this reason the M n cannot be calculated by our 
exact method above. Yet, in certain cases finitely many terms 
of the series are sufficient for a good approximation. 

We first of all express the resonance shift 8h = (h) — co and 
the mean square deviation from the center of the absorption 
peak Ah 2 = (h 2 ) — (h) 2 in terms of the M„, 



8h _ M x 

T = Mb 



Ah 2 _ M 2 M\ 



J 2 



M M 2 



(14) 



There are at least two cases, where these formulae simplify and 
finitely many of the m„ are enough to determine a systematic 
approximation to 8h and Ah. 

The equation for the resonance shift simplifies for small 
anisotropy, |5| <C 1. Since Mo = mo + 0(5), M\ = —mi + 
0(c) 2 ) and, generically, ;«i itself is of order 8 (see ( |T0] >, ( p~3] >) 
we obtain to linear order in 8 



8h 
T 



m\ 
m 



(15) 



In previous work 2 8 the same equation was obtained by a more 
intuitive reasoning. It leads to results which compare rather 
well with experiments. 8 However, some care is necessary with 
the interpretation of (15) . Since m\/8 vanishes at 8 = h = 0, 
it follows that m\ = 8 {ah + b8 + . . . ) with some coefficients 
a,b, whence h/J must be large compared to 8 for ( fT5) to be 
applicable. 

Recall that the higher moments m n , n > 2, are of order 
S 2 . Hence, for the line width there is no simplification for 
small anisotropy, like in ( |15) . But there is another measurable 
quantity which does allow for a systematic small-5 expansion 
to first order, namely the integrated intensity kcoMq, since 



Mq = otq — Jm[ 



(16) 



For the resonance shift it follows to linear order in 8 from |6]) 
and ([15) that 8h(T, co) = -8ca(T,h)\ h=a) . For the line width 
there is no such simple relation between Aco and Ah, not even 
for small 8. 

The representation ( [T3) is a series in ascending powers of 
J/T (with still temperature-dependent coefficients). This can 
be used to evaluate ( fT4| ) asymptotically for high temperatures. 
It turns out that the leading terms in the J/T expansion of 
mi and Jm' 2 cancel each other (8h ~ J^8 — > in the high- 
temperature limit T ^> J) and 



Ah 

T 



V2 
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4T 
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105+9)7 2 +4co 2 
32T 2 
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where CO is the microwave frequency. This formula provides 
a simple means to directly measure the anisotropy parameter 
5. For T — > °° it turns into Eq. (10) of Ref. |9]upon a proper 
identification of parameters. 

In this case as well the momentum-based line width Ah is 
initially slightly increasing from its infinite-temperature limit 
|5|/-\/2 when the temperature is reduced. But our numeri- 
cal data (see crosses in Fig. |4) and the second order term of 
the high-temperature expansion show that Ah behaves non- 
monotonically and decreases again for temperatures lower 
than J. For small temperatures it seems to approach zero lin- 
early. The latter type of behavior is in accordance with field 
theoretical predictions 5 and experimental results!^ We would 
like to stress, however, that there is no contradiction between 
the broadening shown by the upper curve in the first panel of 
Fig. [4] and the narrowing shown by the lower curve. In fact, 
both curves were obtained from the same numerical data set for 
the dynamical susceptibility. The upper curve was calculated 
with |7]), whereas the lower curve was calculated by means of 
( |14) . What is important is that the upper curve can be com- 
pared with our exact results (solid line in the upper panel). The 
good agreement of the crosses with the exact curve creates con- 
fidence in our numerical data. It shows that they are reliable 
when used in integrations. This is a non-trivial statement, since 
our numerical data happen to be noisy and finite-size affected 



at low temperatures (see Sec. V E i. 

We would like to point out that the resonance shift 8 CO /J 
or 8h I J and the line width Aca /J or Ah /J defined in terms of 
moments show a simple scaling behavior. They depend on the 
exchange interaction only through the ratios T/J and h /J. In 
this sense the curves in Figs. T]|4 are universal. 

The method of moments is not only useful for the integrable 
Heisenberg-Ising chain. It may be applied to non-integrable 
spin chains and to two- and three-dimensional models as well. 
The field-dependent moments and the corresponding shifts and 
widths may be accurately calculated by approximate meth- 
ods, since they are determined by static short-range correlation 
functions. For a discussion of the numerical calculation of 
the moments in one dimension see below. In any case, the 
frequency-dependent moments are harder to obtain, since they 
require the calculation of an infinite number of static correla- 
tion functions. 



C. Integrated intensity 

An important quantity in experiments of electron spin reso- 
nance is the integrated intensity. 21 As in the definition of the 
moments we may either integrate over the frequency or over 
the magnetic field. For fixed magnetic field h our definition |5]) 
of the moments mo and nti implies that 



dco°^x'l_{co,h) — n{Jm\ +hmo) . (18) 



As explained in the previous section, in usual ESR experi- 
ments the absorbed intensity is measured as function of h for 
fixed microwave frequency CO. By definition the corresponding 
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FIG. 4. Line widths Aco/J (green) and Ah/ J (red) in the critical 
regime at S = —0. 1 as functions of the temperature. Data from a fully 
numerical calculation for a finite chain Hamiltonian of 24 sites. The 
black cross in the lower panel marks the infinite-temperature limit 
|5|/v2- The solid orange lines are high- temperature expansions of 
Ah I J according to up to first and second order mJ/T. 



In the remainder of this work we shall try to draw at least a 
qualitative picture of how the lines are shaped by considering 
all available limiting cases, where exact results are known, and 
by complementing these with numerical data. In this section 
we review those limiting cases where exact results are known. 
In the following section we describe our numerical calculations. 
Finally, in Sec. [V] we shall consider the line shapes for 8 > 
and T = h = in two-spinon approximation. 



A. Heisenberg limit 

The simplest case where we know the line shape exactly 
is the isotropic case 8 = 5 It may be called the Heisenberg 
limit of the Heisenberg-Ising chain. In this case H is still a 
complicated many body Hamiltonian, but S + commutes with H 
and the time evolution of S + is driven by S z alone (see App.[A|. 



Thus, 5+(f) =e- 1/! 'S+, and 



I(co,h) = 7c5(o)-h)hm(T,h). 



(21) 



This means that there is a single sharp peak, and the absorbed 
intensity is proportional to the magnetic energy hm(T,h) per 
lattice site. This case includes the familiar paramagnetic reso- 
nance (Zeeman effect) for which the magnetization is known 
explicitly, namely m(T,h) = j th ( ^ ) for J — 0. In the general 
case the magnetization must be calculated from solutions of 
linear and non-linear integral equationsPS In our context we 
infer from |9) and ( fTTj ) that m(T,h) = mo = — j<P(o)> i-e. the 
lowest moment mo alone characterizes the line shape. 



B. Ising limit 



integrated intensity is 

/f_V) = j™Jh®x'l_{co,h)=K(oM {T,(o). (19) 

In the paramagnetic regime, where h is large compared to J, 
the integrated intensity is proportional to the magnetization 
m(T,co). The high-temperature expansion of the frequency- 
dependent moment Mq, 

M Q (T,co) = ^ + ^(8-l) + ..., (20) 

following from < fT~3] >, provides another means to determine the 
anisotropy 8 from experimental data. 

III. EXACT LINE SHAPES 

As we have seen the method of moments allows us, at least 
in the field-dependent case, to obtain exact characterizations of 
the resonance shift and the line width for arbitrary temperatures, 
magnetic fields and anisotropy parameters for the infinite chain. 
Unfortunately, it does not teach us much about the actual line 
shapes. Even the most elementary question, how many peaks 
the line comprises, remains generally unanswered. 



The only other limiting case in which the line shape is known 
for all temperatures and magnetic fields is the Ising lirnitPSFor 
the Ising limit we replace the Heisenberg-Ising Hamiltonian 
H — » H/8 and send 8 — » °°. Then 

L 

H^H I =J^s)_ l s). (22) 

;=i 

Due to ([8]l this implies that the moments are replaced as m„ — > 
8~ n m„ for 8 — > oo. 

In the Ising limit the time evolution of S + in (|2]l can be 
calculated explicitly. As we show in App. |A6| this leads to the 
formula 

—X'i-(0),h) = l -{m 2 -m x )8{co-h+J) 

+ (mo-m2)8(co-h) + -(m2+mi)8(a)-h-J) (23 ) 

for the dynamical susceptibility. 

Using this formula we can calculate the moments m„ by 
means of |5]l and verify its consistency. Integrating over (0 we 
obtain indeed mo. The integrals for the higher moments yield 

I mi if n G N is odd 
m„ = < (24) 
W2 if n € N is even , 
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i.e. there are no 'new moments' for n > 2. The three indepen- 
dent moments mo, m\ and m^ correspond to the three 5 -peaks 
in the dynamical susceptibility. They can be calculated by 
means of the 2 x 2 transfer matrix 24 of the Ising chain. Explicit 



T/J = 0.1 



T/J = 5 



expressions are shown in App. A 6 



The Ising limit is important for us, since it is easy to inter- 
pret and since it provides a physical picture for the massive 
phase. The eigenstates of the Ising chain Hamiltonian are ten- 
sor products of local s z eigenstates. We infer from the spectral 
representation, App. |A7| that transitions can occur only be- 
tween states which differ by a single flipped spin. Thus, the 
following transitions in which the chain absorbs the energy AE 
are possible: 



ttt- 

— 1 — 1 — il — 



AE=J-h 
AE=J + h 
AE = h 



The first two correspond to the creation of a pair of domain 
walls (or spinons) in one of the Neel ground states. The third 
one is impossible in the ground states, whence the coefficient 
mo — ni2 in front of the corresponding term in (23) must vanish 
at zero temperature. 

At finite 8 > the operator S~ in the spectral representation 
( |A35| > does no longer induce transitions between eigenstates. 
The three 5 -peaks in (23) broaden which reflects the onset of 
interactions between the spinons. Still, we believe that three 
peaks are characteristic of the massive phase, 8 > 0, at least 
at not too large magnetic fields. This is in accordance with 
our numerical data and with previous numerical work!^ The 
relative height of the three peaks is still approximately well 
described by the relative prefactors of the 5-peaks in (23) . At 
T = in the two-spinon approximation (see below) the central 
peak vanishes due to the same intuitive argument as given 
above. 

Figure[5]shows the dynamical susceptibility (23) in the Ising 
limit as a function of h/J. We visualize the 8 -peaks by con- 
volving them with a Lorentzian of the form 



L{h) 



J 



n h 2 



While in the low-temperature limit (T/J — 0.1) the right peak 
at h = co + J is the highest, in the intermediate-temperature 
regime (T /J = 0.5, 1.0, 5) the relative height of the central 
peak at h — CO increases rapidly with increasing temperatures. 
For even higher temperatures (T/J — 5, 10, 50) the relative 
height of the left peak at CO = h — J increases as well and 
reaches values comparable to those of the right peak. The 
relative heights of the central peak and the right peak are com- 
patible with the numerical data for 5 = 1 shown in Fig. 8 of 
Ref. |6[ The absolute value of the height of the central peak 
differs approximately by a factor of 8, because the authors of 
Ref.|6]plot the function x'L = (X+- + instead of J%+- 

and scale it with a factor of 1 + 8 = 2. Concerning the loca- 
tion of the three peaks as well as the relative and the absolute 
heights of the central peak and the right peak, the curves in 
Fig. 8 of Ref. [6] can be qualitatively explained by the exact 
result (23) for the dynamical susceptibility in the Ising limit. 




FIG. 5. Dynamical susceptibility Jx+- from (23} as function of h/J 
for CO = J. The 5-peaks are convolved with a Lorentzian (25} with 
parameter £ = 0. 127. 



If we use (23} and the formulae ( |A32} for the moments tn\ 
and mi, we can explicitly perform the infinite-temperature limit 
and obtain the function (j)(co — h) defined in the next section 
in (26} and (27) . We find that its 8 -peaks are weighted by 1 /2 
for the central peak and by 1 /4 for the two side-peaks. 



C. High-temperature limit 



In the infinite-temperature limit the dynamical susceptibility 
(25) X+- vanishes identically. This follows, for instance, from 



Eq. ( A37} . From the same equation and from the sum rule 
(A38 1 we obtain the leading high-temperature contribution to 

vC+-> 



con 



x'l_(coJi) = ^$(co-h) + 0(T- 2 ), (26) 



where 



^)(C(J) = — y df e ift "tr{(e ifad "S + )S-}. (27) 

The function is even, non-negative and normalized. Hence, 
it may be interpreted again as a distribution function. 

The function is simpler as compared to X+-- Still, in 
general, we are unable to calculate it exactly. In App. [C] we 
comment on the small-f expansion of the integrand, which we 
have calculated up to the order t 38 , and draw some conclusions. 
In the free Fermion case 8 = — 1 it is possible to calculate it to 
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all orders. The terms sum up to a Gaussian|Sland 

P -(o/J) 2 

0(£») 



5 = -1.0 



(28) 



From this explicit result we can calculate the field- and 
frequency-dependent line widths of the previous section, 



Aco 1 
I 2 



3/2 + 2 {h/J) 4 
(1+2(W) : 



A/i 

T 



1 

V2' 



(29) 



The second equation is in accordance with the high-tempera- 
ture result (17) . Our explicit example clearly shows that %+- 
is asymmetric in h and CO and that the two line widths Aco and 
Ah are rather different quantities. 

We learn from the above discussion that 2T#"_ ((0,h) / KCO 
is a 'good function' . It converges to a normalized function 
which depends only on the difference CO — h for T — > °°. In 
the free Fermion case 8 = — 1 and in the isotropic case this 
function has a single peak at CO = h. From our numerical data 
(see Fig. [6]) we see that this seems to be true for all values 
of 8 between —1 and and even for small positive 8. The 
Gaussian decay for large CO seems to be peculiar of the free 
Fermion point. At all values of 8 which are larger than — 1 our 
logarithmic plots in Fig. [6] indicate an exponential decay. 

When looking for a simple model for such type of line shape 
we found the so-called 'normal-inverse Gaussian', 



N(x|a,j3) 



o/3e^gi(avg+g 



a,/3>0, (30) 



where K\ is a modified Bessel function. It becomes a Gaussian 
in the limit a — > °°, a Lorentzian for a — > and a 5-function 
for j3 — > 0. Its moments can be easily calculated from its 
characteristic function 



>T(jfc|a,j3)= / dcoe lta N(x|a,)3) =e 



For instance, 



a 



(31) 



(32) 



This can be compared with the dimensionless moments of 
the distribution function (j), which follow from |5]) and ( f26) , 



(co 2n )t ,. 4r 



(33) 



Here the right hand side can be easily calculated. If we demand 
that 



(x 2 ) 



CO 



J 2 2 



(x 4 ) 



./ 4 



— - + <5 + <5 2 
2 \2 



(34a) 
(34b) 




FIG. 6. Normal-inverse Gaussian (black lines) as a model for the 
high-temperature line shape, comparison with numerical data (red 
lines). Parameters of the normal-inverse Gaussian as calculated in 
([35j. All panels show J%'J__((0,h) as a function of h/J, left panels 
logarithmic scale, right panels linear scale. Note the different scale on 
the x-axis in the right panel for 5 = —0.1. In general the numerical 
data were obtained for TJJ = 100, CO /J = 0.4, L = 16, and M = 1024 
(see below). For 8 = —0. 1 the chain length was increased to L = 20 
and the resolution to M = 4096. 



we obtain 



a = 



(1 + 5X3-5) 



, P = 



2 V (l + 5)(3-5) 



(35) 



Figu re [6| compares the normal-inverse Gaussian with parame- 
ters (|35| with our numerical high-temperature line shapes. 
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For 8 = — 1 the model line shape is exact. This is no longer 
the case for 8 > — 1 which can be seen by comparing the sixth 
normalized moment of with the sixth moment of N. Still, 
we find it remarkable how well the model line shape fits our 
numerical data. Especially the exponential tails visible in the 
left panels of Fig.[6]have not been fitted to the numerical data. 
The good agreement comes out automatically. On the other 
hand, the deviation of the center of the peak in the right panel 
for 8 = —0.1 in Fig.[6]does not seem to be due to a resolution 
problem of our numerical calculation. We rather attribute it 
to a slight mismatch of the normal inverse Gaussian at small 
anisotropy. 

If we fix the parameters a and /5 of 3^(x\a,j5) according to 
( [35) , we see from ( |34a) that the width, calculated by its second 
moment, behaves as \J (x 2 )^- ~ 1 8 \ . On the other hand, it is 
easy to calculate the peak-to-peak width of the normal-inverse 
Gaussian, which is the distance of its two inflection points. 
Setting the second derivative d 2 N(x\a,fi) to zero and using 
the differential equation defining K\, 



y 2 K'( (y) = (y 1 + l)K l (y)-yK' 1 (y), y = a^x T - 



J3 2 , (36) 



we see that the positive inflection point yo is the solution of the 
equation 



y(y 2 -(al3) 2 ) 
3y 2 -4(aj3) 2 



(37) 



In principle this algebraic equation can be solved numeri- 
cally, but for our purposes the following argument leading 
to an estimation for the inflection point is sufficient. Since 
d y In (Ki (y)) < — 1 for all y > and the left-hand side becomes 
large and positive for y —> 0, °°, the solution yo has to be located 
close to the left of the pole. Hence, we obtain 



< 



(38) 



An upper limit of the peak-to-peak width is therefore given by 



^f=2x <8 2 ] 



(l + 5)(3 + 5) 



(39) 



Accordingly, the normal-inverse Gaussian is an example of 
a distribution for which the width calculated by its second 
moment (~ |5|) and the peak-to-peak width (~ 8 2 ) behave 
asymptotically differently for small anisotropics 8 and can 
therefore differ strongly in value. Fitting experimental data at 
high temperatures to a normal-inverse Gaussian, this offers a 
way to determine 8 independently of the exchange integral J 
from the ratio of the two line widths, 



AppX 



2151 



V^ 2 )^ v/(l + 5)(3 + 5) 



(40) 



If we estimate the peak-to-peak width of Fig. 6 of Ref.l26l 



to Apph s» 2kOe = 0.27 K, we can solve ( |39| ) numerically (J = 
22K) and obtain 5 ~ —0.12. This value is compatible with the 
prediction of Maeda et a —0.15) obtained from a fit for 



the resonance shift based on the data of Ref.[26] The authors of 
Ref.|27] obtain for the same material as in Ref. 26 (LiCuVGj) 
the exchange integral J = 30K. From Fig. 4 therein we can 
read off A pp /i = 1.5 kOe, and obtain together with Eq. ([39) the 
anisotropy 8 « —0.088 which is compatible with the value 
J zz /J « -2K/30K re -0.067 of Ref. El 



IV. NUMERICAL LINE SHAPES 

The numerical approach we are using for the calculation 
of the dynamical susceptibility j;"_(fti,/i), Eq. p), has been 
described in detail elsewhere.^ We will therefore give only 
a short outline of the method and discuss a few special tricks 
beneficial for the present project. 

Starting point of the numerics for finite chains is the spectral 



representation (see also App. A7i 



-E„/T 



-E m /T 



\(m\S-\n)\ 2 8{co-E m +E„) 



LZ 



dys(y + C0,y) (e y > 



y/T . 



(y+(0)/T 



(41) 



which can be written as an integral over thermal weighting 
factors and the temperature-independent function 

s(A:,y) = ^|H^I")| 2 5(^-£m)5(y-£„). (42) 



At first sight the calculation of this function seems to 
require knowledge of all eigenvectors and eigenvalues of 
the Heisenberg-Ising Hamiltonian on a finite lattice. How- 
ever, it is significantly more efficient to rescale all energies 
E — > E = aE + b, such that E g [—1,1], and to expand s(x,y) 
in terms of Chebyshev polynomials of the first kind T(, 



s(x,y) 



M y l mig} {2-8«>){2 - SjoM^Tjjy) 



(43) 



The problem then reduces to the calculation of the expansion 
coefficients which are given by traces, 

IUj = jj dxdys(x,y) T,(x)Tj(y) = rr[S+7K#)S~2>(#)] . 

(44) 

Instead of summing over the whole Hilbert space, these traces 
are well approximated by averages over a few random states. 
The action of T^{H) on an arbitrary state can be quickly eval- 
uated with the recursion relations of the Chebyshev polyno- 
mials. Taking into account symmetries of the Hamiltonian 
(^-conservation and translation), it is thus feasible to calculate 
the flij for systems of up to L = 32 lattice sites and expansion 
orders up to M = 4096 on average hardware. 

Once we have a complete set of expansion coefficients /Jy 
for a given lattice size L, anisotropy 8 and all S z sectors, we 
obtain s(x,y) from Eq. |43jl using fast Fourier methods. In 
Eq. (|4"3"), the damping factors cure the Gibbs oscillations 
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inherent to truncated Chebyshev (and Fourier) expansions, and 
ensure good convergence properties (see Ref.|28 for details). 
Given s(x,y) we can calculate x+- for all temperatures 
T, frequencies CO and magnetic fields h via straightforward 
numerical integration. Changing any of these three parameters 
does not require a new Chebyshev expansion, which is the 
most time-consuming part of the simulations. 

A little more care is required for low temperatures. Here the 
Boltzmann factors put most of the weight on very few states at 
the lower edge of the spectrum. The sums in Eq. ( |4Tj ) should 
then be split up into contributions from these low-energy states 
and from the rest of the spectrum, and the low-energy eigen- 
states should be calculated exactly with Lanczos recursion!^ 
For the data in Figs. [T] and [2] we separated two states per S z 
and momentum sector from the rest of the spectrum. This 
procedure does not increase the overall computation time, but 
the book keeping is slightly more elaborate. 

Another trick improves the precision of the numeric esti- 
mates of the moments m„(T,h) and M„(T,co) of x'L_{(0,h), 
Eqs. (|5]l and ( fT2| ), and of the resonance shifts and line widths. 
Since the expansion coefficients are based on averages over 
random vectors, they are subject to a low level of noise which 
is carried over to ft). Even though the error is hardly 

visible in %'i__{(0,h) itself, it is amplified when %'i__((0,h) is 
multiplied by powers of (a> — h) in the course of the moment 
integration. Large values of (co — h) taken to the power 2 or 
3 then induce noticeable errors in m n (T,h) or M„(T, co). This 
can be avoided by doing the multiplication with (co — h) in the 
space of Chebyshev moments. Consider a one-dimensional 
Chebyshev expansion of a function / : [—1,1] — >M, where the 
expansion coefficients are given by 



form factors, since these are sometimes known exactly. Such 
type of procedure is efficient if sub-classes of form factors can 
be identified which contribute dominantly to the considered 
correlation function. 

For — 1<A:=1+5<1 the Heisenberg-Ising chain is at a 
critical point in the ground state. In this case all form factors 
vanish algebraically in the thermodynamic limit r°" 32 lThus, the 
summation of form factors and the thermodynamic limit do 
not commute. In this case good results for the dynamic struc- 
ture factor were obtained from a summation of form factors 
for finite chains] 33 " 36 ! which requires a considerable amount 
of numerical calculation, though. More recently, an exact 
summation of the leading contribution to the large distance 
asymptotics of two-point functions was obtained in Ref.[37l 

In the massive phase at A > 1 the situation is mathematically 
less involved. Certain classes of multi-spinon form factors 
stay finite in the thermodynamic limit.Ei In calculations of the 
dynamic structure factor,^ it turned out that the two-spinon 
contribution is always dominant at T = h = 0. Here we use the 
results of Refs. [T2"UT4l to discuss the two-spinon contribution 
to the dynamical susceptibility for A > 1. As we shall see, 
the dynamical susceptibility is dominated by the two-spinon 
states only if A is large enough. For A > 3/2 the two-spinon 
contribution I^ 2 \co) to the absorbed intensity amounts to the 
main part of the total intensity I(co,h = 0), but it is marginal 
in the isotropic limit A — > 1 . In the Ising limit A — > °° both 
intensities are identical, I^ 2 \co) — I(co,h — 0), and the result 



of Sec. IIIB is reproduced. 



A comparison of the two-spinon line shapes with numerical 
line shapes calculated for finite chains shows the high quality 
of our numerical data. 



dxf(x)Ti(x). (45) 

i 

Then, the expansion coefficients of xf(x) are 
Tk = J dxxf(x) T,(x) = J dxf(x) T\ (x)Ti(x) 

= \ f_tef(x) (T l+l (x) + Ti_y) = (46) 

Hence, multiplication of the expanded function with the inde- 
pendent variable corresponds to taking a kind of mean value 
in the space of expansion coefficients. The application of 
this procedure to the two-dimensional expansion required for 
(co — h)"x+- ((0,h) leads to a cancellation of noise and to much 
better estimates of the moments m n (T,h), M„(T, co), and of the 
resonance shifts and widths. 



V. TWO-SPINON LINE SHAPES 



Spectral representations such as ( |A35| l, ( |A36[ > have been 
used in the past as a starting point for the approximate cal- 
culation of dynamical correlation functions. For the ground 
state case, Eq. ( |A36| >, approximate calculations can be built 
upon the partial summation of matrix elements of local oper- 
ators between the ground state and excited states, so-called 



A. Line shape 

The dynamical susceptibility x+- decomposes into two 
terms, one for positive and the other one for negative frequen- 
cies. In order to calculate these terms separately we define the 
function 



X((o) 



1 

2L 



d/e i<B, (S+(f)S-) 3 



(47) 



Using the invariance of the Hamiltonian ([TJ under spin flip we 
obtain 



X'i-(co,h = 0) =x(co)-x(-co), 



(48) 



where #(et)) vanishes for CO < and is non-negative for CO > 0. 

The space of excited states decomposes into scattering states 
of an even number of 2n spinons. 38 The precise mathematical 
structure of the space of states of the infinite XXZ chain in the 
massive phase was identified in Refs.fTT1and[39] In order to 
utilize the results of Ref.QTJwe have to adapt our conventions. 
The Hamiltonian 



Hjm 



- J E 



As U s V 



(49) 
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used in Ref . QT| is related to our Hamiltonian ([1} by a uni- 
tary transformation s X j' y i-> (— ) J sJ } and s Z j i-» sj. Under this 
transformation the dynamical susceptibility turns into 



zW = !l /^d/e^OlRV^olO),', (50) 



k=-o 
7=0,1 



where the time evolution in st(t) has to be evaluated by means 
of the Hamiltonian Hjm instead of H, and where the states 
|0)o,i are the two degenerate ground states of Hjm- Here we 
have used 5 ± = Y,k s t' me mvar i ance of the Hamiltonian under 
translations as well as the unitary transformation defined above. 
Now we insert the resolution of the identity into multi-spinon 



states, 



ii 



x \$n---$l)e n ...e l ,jj,e l ...e n ($l---$n\, (51) 



in between the spin operators in ( J50| > and consider only the 
term with n — 2 which is the two-spinon contribution. The 
subindices j refer to the two ground-state sectors and the Q 
are spin indices labeling the scattering states of even numbers 
of spinons. In this language the two-spinon contribution is 



* (2) M = , E I 



jj =0,1 k 

£!,£ 2 =± 



dfe 1 ' 



2 d%i 



f (0| (-) V(0 l&^i)^ J /« (ll , fcko 10)/ • (52) 



Very similar expressions were evaluated elsewhere! 13 l 14 lThe 
time evolution of /(0|(— ) k st(t)\%2,£i) as we U as the remain- 
ing form factors in ( f52] i were obtained in Ref. [TT] Inserting 
those results and calculating the remaining sums and integrals 
we end up with 



X {2 \co) 



k! 0(1 -&)©{& -kf) $l{6) 



(53a) 



Here the factors © in the numerator denote unit-step func- 
tions. The variables 6 and co are related as 



03 A 

0)=-=dn 



2K 







o<e< 



K 



% J 1 ~ ~ 2 
where dn is a Jacobi elliptic function and where 



JK / kK' 

1= — sh 

% \ K 



(53b) 



(53c) 



The anisotropy parameter — (1+8) = (q + q 1 )/2 is related 
to the nome 



q=-exp(-7CK'/K), 



(53d) 



and thus determines the moduli k, k 1 = \J\-k 1 of the elliptic 
integrals K = K(k), K' = K(k'). The remaining theta functions 



in |53) are standard and defined by 



MO) 

y(u) 

(x,y,z) ■■ 



y{q- 2 )y{q- 2 Y 

(q 4 u,q 4 ,q 4 )(u- 1 ,q 4 ,q 4 ) 
(q 6 u,q 4 ,q 4 )(q 2 u- l ,q 4 ,q 4 ) ' 

: n (i-^")- 

m,n=Q 



(53e) 
(53f) 
(53g) 
(53h) 



The line shape of the function x^ is shown in Fig. [7] for 
several values of A. One can observe that the broadened peak 
is very asymmetric for small A. For increasing A the peak 
becomes narrower and more symmetric. The Ising limit A — > °° 
is analyzed in the next section. 



x 




FIG. 7. Two-spinon contribution to the dynamical susceptibil- 
ity Jxl- and corresponding intensity I^ 2 ' as functions of co/JA for 
different anisotropy parameters A>latr = ft = 0. 



B. Ising limit in the two-spinon case 



As the Hamiltonian ( |49] i diverges for A->«iwe rescale all 
energies by the factor A. Replacing, in particular, J by J/A in 
|49|, we obtain the Ising Hamiltonian |22} in the limit A — » 00. 

In the express ion d53"] l for x^ 2 \ the rescaling only pertains 
to the definition ( |53c|i of /, where J must be replaced by J/A. 
All other relations ( |53e| l-( |53h"l i and particularly ( |53b| ) remain 
unaffected, and the Ising limit can be easily performed. Since 
A = (p + p )/2, we conclude that p — > and consequently 
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rescaled / becomes 



l,K-> n/2 as well as K' ->• - \ In (/?) 



7ZA 



sh (nK'/K) 



J p-p- 



2p + p~ 



The 



(54) 



On the one hand the product of the two unit-step functions in 
the numerator of ( |53"j ) ensures that j( 2 ' (a)) = for all CO ^ J, on 
the other hand we show that /* M d(D^f^((D) = n/2 in App. 
Hence, is a 5-function with prefactor 7r/2, 



B 1 



X W{(o) = ^8{(o-J), 



(55) 



which coincides with the result ( |2"3"j ) of Sec. |IIIB| because 
m o = m 2 = for /; = and wtj 1/2 for /; = and T 0. 
For the integrated intensity we easily obtain 



CO 



do) O)X^(co) 



T' 



(56) 



which agrees with Eq. ( p~8^ > of Sec. II C for /i = and m\ = 1/2 



C. Heisenberg limit and integrated intensity 

In the Heisenberg (or isotropic) limit A — > 1 we have p — > 1 
and consequently k — >• 1, — >• 0, /T' — >• 7T/2 as well as K — s- oo. 
Fig. [7] shows that the function %^ tends to zero uniformly. 
This is consistent with the behavior of %+- m the isotropic 
limit, Eq. ( |A24| i, because m(Q, T) = for all temperatures T. 

In order to obtain a measure for the relative contribution of 
to the full susceptibility %+- for all A, especially for the 
isotropic limit A — > 1, we compare the two-spinon contribution 
of the integrated intensity 



/ff(A)= / do)O)^(a)) 



(2), 



with the total integrated intensity 



7 int (A) = / dcocox(co). 
Jo 



We denote their ratio 



r(A) = 



«(A) 

/int(A) 



(57) 



(58) 



(59) 



If we substitute o by by means of ( |53b[ ), the numerator of 
59l becomes 



r ( V(A)- 2k ' KI r /2 de^ 2(0) 
4* (A)-— y d0 ^y 



(60) 



The integral on the right hand side can be easily evaluated 
numerically. Furthermore, we can express the behavior of /> ' 
in the isotropic limit analytically in terms of 1 — p = \/A 2 — 1 — 
(A-l), 



The derivation of this formula and the value of the constant C 
are shown in App. |B 2| 

For the denominator of d59l we use a sum rule and obtain 



2k8J 

/ int (A) = ((sis 



2)0 



(62) 



Using Eqs. (28), (34), and (35) of Ref.@0]the two-point cor- 
relation functions on the right hand side can be expressed by 
integrals which again are easy to compute numerically. We 
obtain 



(■Sl^o- (4*2)0 = ~T + 



dx 



rjch(f) 



3 sin 2 x ch 4 ^ + cos 2 x sh 4 ^ 



j) sin2xshrj 



4(sh 2 f +sin 2 x) 



(63) 



where the parameter rj is defined by A = ch rj . An expansion 
close to the isotropic point rj =0 yields 



/. (A) (7-41n2)K7 4 



(64) 



Accordingly, the total integrated intensity 7j nt (A) tends to zero 
for A — > 1, but nowhere nearly as fast as the two-spinon con- 
tribution /^ t '(A) ~ e -^/ 2 (i-p) \y e conjecture that, close to 
the isotropic point A = 1, all higher spinon contributions X \ 
n>2, are as marginal as J^ 2 ', but in such a way that for the 
full susceptibility % = Ln>i z' 2 "' st iU holds. 

In Fig.[8]the ratio r is plotted as a function of the anisotropy 
A > 1. For A > 3/2 the two-spinon contribution accounts for 
more than 80% of the integrated intensity, for A > 2 for even 
more then 96%. In the limit A — >• °° it rapidly approaches 
100%. Additionally, one can observe in the inset of Fig.[8]the 
over-exponential decay ( |6T] > for A — > 1 . 



< 



1.0 
0.8 
0.6 
0.4 
0.2 
0.0 



FIG. 8. Ratio rasa function of the anisotropy A > 1 . The inset 
shows the behavior of r close to the isotropic point A = 1 , where it 
decays over-exponentially. 
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0.00 
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1.08 



/g (A) -A^ce-W^ (1 + 0(1 -p)). (61) 
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D. Line width and some moments 

In this section we compare the line width of the two-spinon 
contribution to the dynamical susceptibility with the exact line 
width. For this purpose we set T = h = in the expressions for 
cp, CO, and co' in App.[D| Then, all integrals involving functions 
defined as solutions of integral equations vanish and cp, CO, 
and co' are determined solely by their explicit contributions. 
The formula for the line width simplifies to Aco/J = \Jm$fm\ 
which, for T = h = 0, is then expressed in terms of explicit 
integrals by means of Eq. ( fTT| . Equivalently, one may use the 
ground-state results for short-range correlation functions of 
Ref.EU Fig. [9] shows the exact line width as a function of the 
anisotropy for A > — 1 . 



4.5 



3 
<1 



two-spinon line width 
exact line width 
special exact values 



X 




FIG. 9. Two-spinon line width Acq/ J compared with the exact line 
width for T = h = 0. The black crosses marks the special values at 
A= 1,A = 0, andA=-l. 



For the two-spinon contribution the moments mi and m 3 can 
be expressed by the integrals 



(2I)"k'K 







/ dG dn"- 1 


Jo 





(65) 



which can be evaluated to arbitrary precision. The two-spinon 
line width is shown in Fig.|9]as a function of A > 1. As argued 
in the previous section the two-spinon contribution and the full 
susceptibility become identical in the Ising limit A —> °°. An 
expansion in 1 /A of both line widths shows that they agree up 
to the order 1/A 2 . 

When q is a root of unity the moments m\ and m 3 can be 
evaluated!^ Here, we present the results for A = 1, A = 0, and 
A = — 1, respectively, 



m\ (A = 1) 
m 3 (A = 1 ) 



7-41n2 

15 - ' 
2 169 



(66a) 



5 l0 «3)-yC'P)-?CP) 

+ ln(2)(-4-10C(3)+20C(5)), (66b) 



mi(A = 0) 



2(n-2) 



m 3 (A = 0) = 

mi(A = -l) = 1/2 
m 3 (A = -l) =2. 



64-487r-67r 2 + 277T 3 
9tz 4 



(66c) 

(66d) 

(66e) 
(66f) 



The function £ is Riemann's zeta function. The numerical 
values of the line widths at these anisotropics are Aco /J w 
2.00518 for A = 1, Aco/J « 1.84606 for A = 0, and Aco/J = 2 
for A = — 1 . Note that the exact line width as a function of 
A is continuous and non-zero except at the isotropic points 
A = 1,-1 where it is not defined. However, the line width 
can be continued continuously at these points which yields the 
curve plotted in Fig. [9] 

E. Comparison with numerical calculations for finite chains 

We can now compare the two-spinon contribution to the 
dynamical susceptibility % ■< Eq. ( |53] l, with the full suscepti- 
bility obtained numerically by the method described in Sec. IV 
For this purpose we use numerical data for x+- as a function 
of CO for chains of lengths L = 24 and L — 32 at h — and 
A = 2. We shall indicate the length by a s ubscr ipt and briefly 
write x'l- Based on the discussion of Sec. V C we expect that 
the two-spinon contribution amounts to the main part (~ 96%) 
of the full susceptibility x+- m me limit T — > 0. For larger 
temperatures we expect deviations. 

We have to comment on the limit T — > 0. For finite chains 
the ground state with energy £0 is non-degenerate and for 
L mod 4 = carries momentum q — 0. Yet, the gap AE = 
E K — Eq to the low -lying q = 7T-state is very small compared to 
the gaps to all other states. In the thermodynamic limit these 
two states degenerate and stay separated from the rest of the 
spectrum. Hence, for a better comparison with the two-spinon 
contribution, we averaged over the q — 0- and q = 7T-states in 



our numerical calculation shown in Fig. 10 

As one can see, at low temperature the dynamical suscep- 
tibilities X24 anc l X32 ( re d li nes ) consist of a multitude of 
narrow peaks. This peak-structure is due to the finiteness 
of the chain and the small number of eigenstates that con- 
tribute to the T — response. The Chebyshev expansion 
approach, whose resolution is inversely proportional to the 
expansion order M, can then distinguish all contributing matrix 
elements. With increasing temperature the Boltzmann factors 
exp(—E„/T) — exp(— E m /T), see App. A7 suppress fewer 



states, and eventually the density of the peaks becomes higher 
than the numeric resolution. The dynamical susceptibilities 



then evolve into smooth curves, as is illustrated in Fig. 1 1 and 
also by the high-temperature data in Fig. [6] 

Similar behavior occurs for increasing lattice size L. Com- 



paring the two panels of Fig. 10 one observes that the peak 



structure of x'l becomes tighter for larger L. Additionally, the 
heights of the peaks decrease. Although x'l ( re d nnes ) an d X^ 
(dark-blue line) do not look alike, the integrals of these two 
functions (orange and light-blue line) match very well. This 
indicates that the step function J a dot)' %'l(G)') (orange line) 
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JxmH at T I J = 001 

integrated 
integrated 



uj/J 



FIG. 10. Two-spinon contribution compared with the dynami- 
cal susceptibility Jx'l of finite chains as functions of CO /J for A = 2 
at T = h — 0. L — 24 in the upper panel, L = 32 in the lower panel. 



converges uniformly to f® dco' x+-((o')- For the thermody- 
namic limit L-)oowe expect that the peak-structure smears 
out and the dynamical susceptibility becomes a smooth curve 
akin to the two-spinon contribution . 

We conclude that at low temperature the numerically calcu- 
lated line shape of the dynamical susceptibility shows strong 
dependence on the size of the chain and on its finite-size spec- 
trum. By way of contrast, the moments of the dynamical 
susceptibility and other integrals over the whole range of fre- 
quencies seem to be well approximated by our numerical data. 



VI. CONCLUSIONS 

The Heisenberg-Ising chain considered in this work is a 
prototypical model of a quasi one-dimensional anisotropic an- 
tiferromagnet. Its collective spinon excitations are created in 
pairs. Their spectrum is a scattering continuum characteristic 
of one-dimensional interacting systems. In microwave absorp- 
tion it becomes visible in a broadening of resonances away 
from the isotropic point. 

Within the linear response theory the absorbed intensity is 
basically equal to the imaginary part of the dynamical suscep- 
tibility multiplied by the absorption frequency. Although the 
model is exactly solvable as long as the magnetic field is di- 
rected along the axis of magnetic anisotropy, the calculation of 
such type of dynamical correlation functions at finite fields and 




uj/.J 



FIG. 11. Two-spinon contribution Jx^ (dark-blue line) compared 
with the dynamical susceptibility JX24 ( re d li nes ) as functions of 
co I J for A = 2, h = and L = 24 in all panels. The temperature is 
increased from T/J = 0.01 to T/J = 2.5. Note the different scales at 
the vertical axis. 



temperatures is still beyond the possibilities of contemporary 
theoretical methods. Due to recent progress in the calculation 
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of static short-range correlation functions, however, it became 
possible to calculate certain global characteristics of the spec- 
trum: the field-dependent moments of Sec.[ll]that determine the 
average absorption frequency (the resonance shift 8 CO) as well 
as a field-dependent line width Aco at arbitrary temperatures 
and magnetic fields. 

In this work we compared the exact data for the resonance 
shift 8 co and the line width Aco with data extracted from a nu- 
merical calculation of the dynamical susceptibility for chains 
of finite length. We used the exact data to improve the numer- 
ical calculation and to validate the quality of the numerical 
data. When looking at the numerical data for the susceptibil- 
ity they appear spiky and finite-size dependent. In any case, 
they look rather different from the smooth curves obtained for 
the two-spinon contribution to the dynamical susceptibility of 
the infinite chain at T = h = 0, which is quasi-exact at large 
enough anisotropy. However, and this is an important part 
of the moral of our work, the picture changes if we look at 
integrated quantities. The integrated susceptibility for L = 32 



in Fig. 10 appears already rather similar to the integrated two- 
spinon susceptibility. When turning to moments the picture 
becomes even better. The numerical finite-chain data for the 
width Aco in Fig.[2](where they are compared with the exact 
data for the infinite chain) look as if they are almost free of 
finite-size corrections. This gives us confidence that our data 
for the frequency-dependent line width Ah, calculated from the 
same numerical data set for the dynamical susceptibility, are 
reliable as well. 

It was a big surprise for us that Ah shows the opposite tem- 
perature dependence in the critical phase as A© (see Fig.|4|. 
The line width Aco increases as the temperature decreases, but 
Ah decreases. This markedly distinct behavior can be attributed 
to the asymmetry of the dynamical susceptibility in CO and h. 
Not only the temperature dependence of the two measures of 
the line width defined by the two types of moments m„ and 
M n is different, but also their absolute values. We observe that 
Ah < Aco. 

In general, the peak-to-peak width, usually measured in 
experiments 3 4 and decreasing with temperature, cannot be 
extracted from our numerical data, since they are not smooth 
enough at low temperatures. At high temperatures, however, 
where we can use the normal-inverse Gaussian as a model 
for x+- (<Dj A) I CO, we find a peak-to-peak width A/! pp which is 
again smaller, A/! pp < Ah < Aco, and whose magnitude seems 
to be compatible with experiments. 

The conclusion for theoretical attempts to extract the line 
width from approximations to the dynamical susceptibility is 
that the seemingly simple and intuitive concept of a line width 
is rather delicate. The peak-to-peak width, popular in the anal- 
ysis of experimental data, is shape dependent and influenced 
by a priori assumptions on the line width. By contrast, our 
moment-based line width Aco is not based on assumptions 
about the shape of the spectral lines and can be calculated ex- 
actly for the Heisenberg-Ising chain. It is, moreover, universal 
with respect to a scaling of all quantities with the exchange 
interaction J. For these reasons we are curios if it will be 
possible in practice to obtain Aco from experimental data. 

This will depend on how well background and noise can 



be separated from the signal. From clean data one could even 
directly extract the moments m\, ma, m?,, . . ., defined in ((5), 
which would mean to directly measure certain short-range 
correlation functions ranging over 2, 3, 4, . . . lattice sites. 

As opposed to the line width the resonance shift is expected 
to be a more robust quantity. We expect our results for 8 CO to 
compare rather directly with experimental data as long as the 
observed line shapes are not too much asymmetric. In the latter 
case the definition |6|) should be taken seriously and should 
be used to calculate the average absorption frequency from the 
experimental data. 

From the two-spinon result for the absorbed intensity (see 
Fig. [7j it can be seen that the spectral lines at low temperatures 
can be expected to be broad and asymmetric. In principle, the 
amount of asymmetry of the lines can be expressed in terms 
of the higher moments ni4, ms of the dynamical susceptibility. 
And, in principle, these higher moments can be calculated 
exactly at any temperature and magnetic field, which we leave 
as project for future research. Another interesting project for 
the future may be the calculation of the T — dynamical 
susceptibility in the critical regime by means of form factors 
in the finite volume, in analogy with the work of Refs. [3314361 
on the dynamical structure factor. 
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Appendix A: Absorption of energy in quantum spin chains 

In order to keep this work self-contained we include a sum- 
mary of the linear response theory of energy absorption and its 
application to quantum spin chains. 



1. Time evolution of the statistical operator 

We consider a quantum system with Hamiltonian J{ pos- 
sessing a discrete spectrum (E n )™_ and corresponding eigen- 
states {|n)}" = o- At time to we adiabatically switch on a time- 
dependent perturbation V{t). We are interested in the time 
evolution of the system, assuming it was initially, at times 
t < to, in an equilibrium state described by the statistical oper- 
ator 



Po = ^£e r\ n )(n\ 



(Al) 



of the canonical ensemble. We denote the temperature by T 
and the canonical partition function by Z. 
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Let U (t) the time evolution operator of the perturbed system, 

idtU{t)=(3i + V{tj)U(t), £/(f )=id. (A2) 

Under the influence of the perturbation the state \n) evolves into 
\n,t) = U(t)\n), and the statistical operator at time t becomes 

P( f ) = \ 5>~^MM = U(t)p U- l (t). (A3) 
We define 

/?(?)= e iM '(p(0-po)e _i:K ', (A4a) 
W(t) = e m 'V{t)e- mt . (A4b) 

Then 

i#J?(f) = He i ^C/(f)p (e i:K 'f/(f)) _1 

= [W(t), e m 'p(t) e - mt ] = [W(t),R(t)+p ] . (A5) 
Since R(to) = by construction, we obtain 

R(t) = -i fdt'\W(t'),R(t')+po]. (A6) 
J t 

This Volterra equation is an appropriate starting point for a 
perturbation theory. Assuming that W(t) be small we conclude 
that 

R(t) = -i f dt'[W(t / ),p } + O(W 2 ), (A7) 
i.e. to lowest order in W 

p(t)=p -ic- mt fdt'[W(t'), Po }e i:Kt . (A8) 

J to 

This is the statistical operator in Born approximation. In the 
following to will be sent to — °°. 

2. Time evolution of expectation values 

Using ( | A8] > we can calculate the time evolution of the expec- 
tation value of an operator A due to the perturbation. Writing 
A{t) = e'^'Ae -1 "^ and using the invariance of the trace under 
cyclic permutations we obtain 

8(A) T =t r {(p(t)-p )A} 

= -if dt'\i{[W(t'),p }e m 'Ae- mt } 

= -i f Jt'([A(t-t') : V(t')]) T . (A9) 

Here = tr{po } denotes the thermal average. A typical 
example of a perturbation, which will be relevant for our discus- 
sion below, is a classical time-dependent field h a (t ) coupling 
linearly to operators X a , 

V(t)=h a (t)X a 1 (A10) 
S(A) T = -i j^dt'([A(t-t'),X a ]) T h a (t'). (All) 



3. Absorption of energy 

The absorbed energy per unit time is 

d ^ = ±tx{( P (t)- Po )(X + V(t))} 

= -itr{pi + V(t),p(t)]CK + V(t))} 
+ tx{(p(t)-p Q )V(t)} 

= tr{(p(0 -po)V(t)} = 8(V(t)) T . (A12) 

Here we used (jA2j, ( |A3[ > in the second equation and the cyclic 
invariance of the trace in the third equation. Assuming that 
V(t) is of the form ( |A10| > and using ( |A1 1 1 > we obtain 

^=-1^ At' {[X a {t-t'),xP]) T h a (t)hP{t'). (A13) 

4. Application to quantum spin chains 

Let us now apply the above formalism to the Hamiltonian 
of the Heisenberg-Ising spin chain in a longitudinal static mag- 
netic field of strength h, 

^ = J t (4-1^ + ^-1^ + (1 + - hS " ■ ( A14 > 

7=1 

We perturb the spin chain by a circularly polarized electro- 
magnetic wave propagating in z-direction. We assume that the 
wave length is large compared to the length of the spin chain 
and idealize this assumption by setting the wave number k = 0. 
Then the magnetic field component of the wave is 

(cos(cof) \ 
-sin(fi)f) , A >0. (A15) 
/ 

It couples to the total spin as 

V(t)=h a (t)S a . (A16) 

Thus, 

f = -if_^'([S a (t-t'),S%h a (t)hP(t') 

= ^V{e^-0<[ 5 +(0,S+]) r 

^e-^ 2 '- f \[S-(t') 1 S-}) T +^''([S + (t'),S-]) T 

- e - i(0t '([S-(t'),S + ]) T }. (A17) 

The ability to absorb radiation is a material property. Hence, 
we generally expect the absorbed energy per unit time to be 
proportional to the number of constituents of a physical system 
and to diverge in the thermodynamic limit. In order to define 
a quantity that truly characterizes the material and is finite in 
the thermodynamic limit we should therefore normalize by the 
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average intensity A 2 of the incident wave and by the number 
of lattice sites L. Further averaging the normalized absorption 
rate over a half-period Tt/(0 of the applied field, we obtain the 
normalized absorbed intensity 

CO fa . dE 
dt 



CO 

ALJ- 



dte m {[s + {t),s-]) 1 



(A18) 



Introducing the familiar (imaginary part of) the dynamical 
susceptibility per spin, 

x'l-(co,h) = ^ L fjte"»'([S + (t),S-]) T , (A19) 



the normalized absorbed intensity reads 



CO 



(A20) 



which is Eq. ([3) of the main text. 



5. The isotropic chain 

The full dynamical susceptibility can only be calculated 
in certain special cases. In order to understand its behavior 
qualitatively we first of all consider the isotropic point 5 = 
of the Heisenberg-Ising chain. In this case 



pK,S] = -h[S z ,S], 



(A21) 



and the Heisenberg equation of motion for the total spin S can 
be solved, 



= ipK.S*] = -iA[S !Z ,5 ± ] = . 
^S ± (t)=e* ih, S ± , 
and S z (t) — S z . Hence, the total spin behaves as 

s(o 



cos(ht) sin(/zf) 
-sin(fe) cos(fe) | S. 



(A22) 



(A23) 



It rotates clockwise about the z axis. 

On the other hand, inserting ( |A22| > into ( |A19| > we obtain 

Z+_(<M) ! ' 



dte ii( °- h >([S + ,S-}) T 



2L J- 
2n8{co-h)m{T,h) 



where m(T,h) = (S z )t/L is the magnetization per lattice site. 
The corresponding normalized absorbed intensity is Eq. pTj ) 
in the main body of this article. 

Comparing ( |A15| l, ( |A23[ > and ( |A24[ ) we interpret the absorp- 
tion of energy as a resonance between the rotating field of the 
incident wave and the rotating total spin of the chain, both 
spinning clockwise with frequency CO = h. If we are off the 
isotropic point 8 = of the Hamiltonian ( |A14| > we may expect 
that energy is transferred from the 'coherent motion of the 
total spin' to 'other modes', causing a damping of the spin 
precession and hence a shift and a broadening of the spectral 
line. 



Remark. In other treatments of the same problem the incident 
wave was considered to be linearly polarized, leading to an 
absorbed intensity 

co izh 

^x"Aa,h) = —m{T,h)(8{co-h) + 8{co + h)). (A25) 

This can be understood by taking into account that a linearly 
polarized wave can be decomposed into a superposition of two 
circularly polarized waves of opposite circular polarization. 
For this reason the two spectral lines in ( |A25[ > are, in fact, 
one and the same, when either the right circularly polarized 
wave has frequency CO or the left circularly polarized wave has 
frequency — co. 



6. The Ising chain 



For the Hamiltonian Hj of Eq. p2| ) the time evolution 
e 1 ' adw '5 + can be calculated explicitly. One easily proves by 
induction that 

ad^S+=y"£(^ 1+ ^ +i r4 (A26) 

for all non-negative integers n. Furthermore 

{s z _, +s z , , -, if n is odd, 
2 + 2Sj_ j s j + j if n is even 



for all n € N. It follows that 



e ifad "/S+=S + -A- 



-Ut 



(A28) 



where 



A = t ( j + 2*5-1*5+1 )4 . B = t (4- 1 + s 5+i )4 ■ 

(A29) 

Inserting ( |A28| > into the definition of the dynamical suscepti- 
bility (|2| in the Ising limit we obtain 



(A24) Z+-(<M) = ^{([S + -A,S-]) T 8(co-h) 



+ -([A+B,S-]) T 8(co-h+J) 

+ 1 -([A-B,S-]) T 8(co-h-J)Y (A30) 

The coefficients in front of the 5-functions can be easily ex- 
pressed in terms of the moments in the Ising limit, 

^([A,S-]) T =m 2 = l -(s\ +4*f44) r , (A31a) 



1 

1L 



([B,S ]) T = -mi = 2(s\s z 2 )-i 



(A3 lb) 
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For these correlation functions explicit expressions in terms of 
h and T can be obtained by means of the 2x2 transfer matrix 
of the Ising chainj^l 



<4> 



where 



2 A /sh- 



2( h_\ 
,2T/ 



>i/r 



(s 

4 v 



e J/T 
2 



i/r 



\2r' r 



ch(x) - Wsh z (x) +e- v 



ch(x) + Wsh z (x)+e- v 



and 



(444) T = (4)r(2(44)T - (44 



(A32a) 
(A32b) 

(A33) 
(A34) 



Inserting ( |A3 1 ] > into ( |A30] > we obtain Eq. (J23j of the main 
text. 



7. Spectral representation of the dynamical susceptibility 



The dynamical susceptibility has the spectral representation 



x\(m\S-\n)\ 2 8(eo-E m +E„) (A35) 



following from ( | A 1 9| >. Here the E n are eigenvalues of J£, 
i.e. they include the dependence of the magnetic field . 

Let CO > 0. Then the only non-zero terms under the sum are 
those with E m > E„ and those are positive. Hence, j"_(o)) is 
non-negative for CO > 0. Similarly, %'!_ (to) is non-positive for 
© < 0. It follows that I(co) is non-negative as was expected 
on physical grounds. The spectral representation simplifies for 
T -> 0, 

^_( <0 , / , ) ^_^_£{| (n |5-|g)|2 5 ( fi) _ J E„ + ^) 



-|te|S»| 2 S(©-£ g +£ n )}. (A36) 

Here we average over the n g degenerate ground states \g). In 
general, %+- ls non- vanishing even for vanishing magnetic 
field. 



8. Alternative representation and a sum rule 

From the spectral representation ( |A35[ ) we infer that 



K 



LZ 



l-e-T)Y,e-^\(m\S-\n)\ 2 8((D-E,„+E„ 



2L 



l-e-r) / dts im, (S + (t)S-) T . (A37) 



This representation immediately implies the sum rule 

dto X+-(co,h) 1 l^ 1 + _ 

,2^T^F = 2L (5 S > r = 2^ 1>r " 

(A38) 

Since the correlation functions (s^sJ +1 )t decay exponentially 
in j in the thermodynamic limit, the integral on the left hand 
side exists for L — > °°. If %'l((0,h) is analytic in the ther- 
modynamic limit, it follows that the constant term in the 
Taylor expansion at CO — must vanish. Then the intensity 
l((0,h) = CO%'l_(cQ,h)/2 has a double zero at co = 0. Hence, 
it must have at least two maxima, one for CO > and another 
one for CO < 0. 



Appendix B: Technical details of the two-spinon calculations 

1. Integrated susceptibility in the Ising limit 

We calculate the integral of the susceptibility %^ over posi- 
tive frequencies CO in the Ising limit A — » «>, i.e. p — > 0, 

Jo n Jo dn(^e) # 2 (0) 



where we have inserted q33p and substituted CO by via 
the relation ( |53b| i. Due to K —t n/2, k' -» 1 we obtain 
dn(2#0/7z;) -> dn(0) ->■ 1, d„(fl) 1 for all e (0,71/2) 
and thereby 

^ dcox [i] {co)^ \ r /2 d9^(9). (B2) 
Jo 2 Jo 

The limit p — ^ for the function #a is more complicated. Using 
the definitions |53f]l-(|53h} of #a and the relation 



in( n (i-zpi-ptr)) 

\n t ,...,n,„>0 ) 



y i i. 



, (B3) 



we obtain 

y(z 



Y(q~ 2 ) ' r \ sh(2«£)ch(«£) n / ' 

where z — e 2 ' e and e = nK' /K, and therefore 



^ 2 (0)=exp -£ 



ch(2«e)cos(4n0)-l e" e 
! sh(2«e)ch(«e) n 

We now convert dFJSb into the form 



(B5) 



00 ^ 2/7 

ln(# A (0)) = In [2sin(20)] + £ - -f-^ sin 2 (2«0) 



1 p 2n 



I ^2n 



- 2 .?.»ti 



(1 + P 2 ") 2 



cos z (2«0) (B6) 
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and see that ~&a(9) 2 sin (20). We finally obtain 



Appendix C: Short time expansion at infinite temperature 



d(ox {2) {(») 2 ( n 2 AO sin 2 (20) = - . 

o Jo 2 



(B7) 



2. Integrated intensity and Heisenberg limit 

We analyze the behavior of the two-spinon contribution of 
the integrated intensity d60|, 



$(A) = 



2k! KI [*/ 2 , a 0?(0 



d0 



K Jo 02(0) ' 



(B8) 



in the isotropic limit A— > 1, i.e. p — > 1. On the one hand, it is 
known that i?„(0) ch (2K/n0) -)• », yt' -> and / -> 0, on 
the other hand we have K ^ °° and i^(0,/?) — > oo. In order to 
obtain more explicit results we define the function 



/f^0) =lim/tf A 2 (0). 



which can be calculated as 



-12f'(-l) 



2 4/3 



-jc sh (jc) 



-2ix/n 



x exp i—f dx' In 

[ ft Jo 



where y/ is defined by 



y/(x) = 



r(i -k)r(i/2+k) 
r(i+k)r(i/2-k) 



Hence, we obtain 



p-h \ k' J p-n \ n J Jo j?„ 2 (0) 



A straightforward calculation yields 



(B9) 



(BIO) 



(Bll) 



, =2,6471... (B12) 
o ch jc 



In the limit of infinite temperature the dynamical correlation 
functions reduce to traces of the considered time-dependent 
operators, and it is tempting to express the time evolution in 
terms of nested commutators, evaluate the traces, and obtain 
the leading terms of a Taylor series near time t = 0. For the 
spin-pair correlation function (s z n (t )sq(0)) such series have 
been calculated to orders as high as ? 30 already two decades 
agoF^Unfortunately, these series converge slowly and with the 
accessible number of terms precise values of the correlation 
functions can only be obtained for rather short times t < 5. 
Nevertheless, we performed such an expansion for tr[S + (f)S~], 
which is at the core of the function <p(co) defined in Eq. |27li. 
Using the highly efficient computer algebra program FORM, 
we computed the series up to the order f 38 . The first few terms 
read 



1 



2,2 



2Tr[S+(r)S-] wl--(A-lp 



2\A 



96 



(A-l) 2 (3-2A + 2A 2 )r 



1 



4\,6 



11520 



(A- 1) (30— 15A + 50A 2 - 16A 3 + 8A 4 )f 



(CI) 



r4-A 

224T f 



Taking the logarithm we find 

ln(2Tr[S + (f)S~]) « -~(A- 1)V +A(A- l) 2 

15-140A + 76A 2 -8A 3 fi 

H 1 

+ 2 4 6! 

56 - 2730A + 10948 A 2 - 8792A 3 + 2256A 4 - 136A 5 g 

H i f 

2 6 8! 

(C2) 



and we can immediately read off the two known results for 
A = and A = 1, — f 2 /4 and 0, respectively. From the terms in 
the square bracket we were only able to sum the A-free contri- 
butions and the terms of highest order in A. The coefficients of 
the A-free terms read 4, 15, 56, 210, 792, . . . , and correspond 
to 



2k 

1-1, 



with k = 2,3,... leading to the Bessel function 1%. The 
coefficients of the terms with the highest power of A for each 
power of t read 1, 8, 136, 3968, 176896, . . . , and constitute the 
expansion coefficients of tan(;t:) 2 , 



In (k>) = 41n [(p,p 2 )] -41n [{- PlP 2 



p-y\ 7t 



1; . +21n2+ — + Q(l-p). (B13) 



tan(x) 2 



1 



136 



3968 



2 2! 4! 6! 8! 
Summing these two sub-series of (|C2[), we obtain 



(C3) 



(2) 

The asymptotic behavior of /> / therefore reads 



I®(A)-^Ce W-p) (l + 0(l-p)), (B14) 



where C = Aq^^k w 124,6. 



ln(2Tr[S + (f)S~]) w 4A(A- l) 2 / 2 (f) - (A- l) 2 (A/2)f 2 
2(A-l) 2 ln[cos(fA/2) 



A 2 



(C4) 



This function does not have much value as an approximation 
of the considered infinite-temperature correlation function, but 
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it contains at least all known exact results, 



r exp(-f 2 / 4 ) fbrA = 0, 
Tr[S+(f)S-] -> { i forA=l, 
i(l+cos(f)) forA^oc 



t = f A finite. 

(C5) 

The last l ine corresponds to the Ising limit discussed in 
and the Fourier transform of 1 +cos(f) gives the 



Sec. 



11IB 



three 8 -peaks of (j), Eq. fl27) , where in the limit T —> °° the 
side-bands have half the weight of the central peak. 



terms of solutions of non-linear and linear integral equations. 
They were termed 18 the physical part of the problem, since the 
physical parameters like temperature or magnetic field enter 
solely through these functions. We shall provide their defini- 
tion only for the massless case — 1 < 8 < O.^The definitions 
for the massive case can be found in Ref. 16 



Appendix D: Physical part of the static correlation functions 

The functions <p, ft), and ft)' that determine all static corre- 
lation functions of the Heisenberg-Ising chain are defined in 



First of all let us define a basic pair of auxiliary functions as 
the solution of the non-linear integral equations 



lnb(x) 
lnb(x) 



7th 2 ft/ sin (y) 

'2(n -y)T ~ /ych(ftx/y) 
nh 2 ft/ sin (y) 



2(ft-y)/ /ych(ftx/y) 



with kernel 



F(x) 



dk- 



dy 

-=o 2ft 



F(x-y)ln(l + b(y))- 



dy 
2ft 



F(x-v)ln(l+b(v)) 



sh((f -y)k)e ikx 



2sh((ft-y)§)ch( 



(D2) 



Here we have introduced a parameter y which provides yet 
another parameterization of the anisotropy, 8 = cos(y) — 1. 
Eq. ( |D1| > is valid for < y < ft/2 meaning that — 1 < 8 < 
0. Below we shall also use rj = iy. Note that the physical 



dy 

_oo2ft 

dy 
2ft 



F(x-y + f}-)ln(l + b(y)), (Dla) 
F(x-y-r]-)ln(l + b(y)) (Dlb) 



r 



parameters temperature /, magnetic field h, and coupling / 
enter only through the driving terms of Eqs. ( |D1[ ) into our 
formulae. 

Except for the auxiliary functions b and b we need two more 
pairs of functions g^ and g'^ in order to define <p, ft), and 
ft)'. Both pairs satisfy linear integral equations involving b 
and b, 



and 



g'f ) (x) = i f sech 



r 



dy F(x-y) 
2ft l + b-'Cv) 



— 1 , 



-oo2ft 



i + b" (y) 



2n l + b '(y) 

dy F(x — y — tj~) 
-„o2ft l + b-\y) 



-M)-f 



l ) W=(f(x- M )- 



^)sech(^ 

r dy F{x-, 
J-oo 2ft l + b- 1 

-§)sech(^%^ 



r~ dy D{x-y) j +) r dy Djx-y + r,-) ( ) 

r 7- o2ftl + b- 1 (y)^ [y) r 7-»2ft i + b-\y) M 

Hy)^ w Lift 1+F i w W ' 



dy Ffr-y) (_) /™ dy F^-y-T?-) (+) 

2ft 1 ^pr 1 AA g " W A- 2ft l + b-Hv) W ' 



(D3a) 



(D3b) 



(D4a) 



(D4b) 
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where 



of the XXZ chain. The function o'(jUi , jUa) is defined as 



D(x) 



sin(fcc)sh(f )ch((f-y)Jt) 
4sh 2 ((^- 7 )|)ch 2 (f) 



(D5) 



The functions co(ji\fl2), (*>' '(flu M2) and <p(ju) that deter- 
mine the explicit form of the correlation functions of the XXZ 

chain can be written as integrals involving b, b, g^ and g'^. 
The function 

<Km) : 



cbc 



2(w-y)i 



l + b '(x) 



(D6) 



determines the magnetization m(T,h) — —49(0) which is the 
only independent one-point function of the XXZ chain. The 
function 



£»(Mi>M2) 



dfc 
dx 



sh((B: -7)§)cos(fc(/Xi -^ 2 )) 
i*(f)ch(¥) 



ych 
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+ b l ( x ) l + b \x) 



with 



K(fl) = cth(jU - 77) - cth(;U + 77) 
also determines the energy per lattice site 

(^+^_ 1 ^ + A^._ 1 ^.) r = sh(T ? )fi)(0 ) 0)/4 



(D7) 



(D8) 



(D9) 



0)'(Ml,Al2) = §^ (+ Vl-M2) 



dk 



dx 

— ych^-^ 



dx (x-jU 2 ) 
7 ch(^ 



2ith(f)ch 2 (f) 
l + fa-'W 1+F'(x) 

rit } W si? to 



l + b l (x) i + b \ x ) 



(D10) 



where 



= cth(ji-T/) + cth(/i+T/), 



(Dlla) 
(Dllb) 



For the calculation of the moments in Sec.|Il]the non-linear 
integral equations for b and b as well as their linear coun- 
terparts for g^ and g'j? were solved iteratively in Fourier 
space utilizing the fast Fourier transformation algorithm. The 
derivatives of gj, and g'j? with respect to /x, needed in the 
computation of the respective derivatives of <p, co, and co' sat- 
isfy linear integral equations as well, which were obtained as 

derivatives of the equations for and g 1 ^. Taking into 
account derivatives is particularly simple in Fourier space. 
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